#######################################################################################
# Summary of Simulation Results: Monte Carlo Study 2 
#################################################################################

# Function to Average Across S=1000 Simulations

MeanOutFun <- function(IN_FILE, OUT_FILE){

result.array <- readRDS(paste(IN_FILE, ".rds", sep=""))
mean.array <- matrix(NA, 6, 8)
rownames(mean.array) <- rownames(result.array)
colnames(mean.array) <- c("NumRob", "TP", "FP", "TN", "FN", "CSign", "RMSE", "Corrs")
for(i in 1:6){
  for(j in 1:8){
    
    mean.array[i, j] <- mean(result.array[i, j, ], na.rm=T) 
    
  }
}

#write.table(round(mean.array, 2), file=paste("/Users/richardtraunmuller/Dropbox/PSRM_Replication/", OUT_FILE,  ".csv", sep=""), sep=",")
print(mean.array)
}

##################################################################################################################################
# Table 2: Probability of False Positive and False Negatives When Variables are Highly Correlated

p <- 1

n1 <- round(1-MeanOutFun("one_runif_result_mat_1", "one_runif_means_25")[c(1:2, 6, 3:5), c(2)]/p, 3)
n2 <- round(1-MeanOutFun("one_runif_result_mat_2", "one_runif_means_50")[c(1:2, 6, 3:5), c(2)]/p, 3)
n3 <- round(1-MeanOutFun("one_runif_result_mat_3", "one_runif_means_75")[c(1:2, 6, 3:5), c(2)]/p, 3)
n4 <- round(1-MeanOutFun("one_runif_result_mat_4", "one_runif_means_10")[c(1:2, 6, 3:5), c(2)]/p, 3)
n5 <- round(1-MeanOutFun("one_runif_result_mat_5", "one_runif_means_125")[c(1:2, 6, 3:5), c(2)]/p, 3)

p1 <- round(1-MeanOutFun("one_runif_result_mat_1", "one_runif_means_25")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p2 <- round(1-MeanOutFun("one_runif_result_mat_2", "one_runif_means_50")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p3 <- round(1-MeanOutFun("one_runif_result_mat_3", "one_runif_means_75")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p4 <- round(1-MeanOutFun("one_runif_result_mat_4", "one_runif_means_10")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p5 <- round(1-MeanOutFun("one_runif_result_mat_5", "one_runif_means_125")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)

n <- t(rbind(n1, n2, n3, n4, n5))
p <- t(rbind(p1, p2, p3, p4, p5))

tab.2.a <- rbind(n[1,], p[1,], n[2,], p[2,], n[3,], p[3,], n[4,], p[4,], n[5,], p[5,], n[6,], p[6,])

#
p <- 5

n1 <- round(1-MeanOutFun("five_runif_result_mat_1", "five_runif_means_25")[c(1:2, 6, 3:5), c(2)]/p, 3)
n2 <- round(1-MeanOutFun("five_runif_result_mat_2", "five_runif_means_50")[c(1:2, 6, 3:5), c(2)]/p, 3)
n3 <- round(1-MeanOutFun("five_runif_result_mat_3", "five_runif_means_75")[c(1:2, 6, 3:5), c(2)]/p, 3)
n4 <- round(1-MeanOutFun("five_runif_result_mat_4", "five_runif_means_10")[c(1:2, 6, 3:5), c(2)]/p, 3)
n5 <- round(1-MeanOutFun("five_runif_result_mat_5", "five_runif_means_125")[c(1:2, 6, 3:5), c(2)]/p, 3)

p1 <- round(1-MeanOutFun("five_runif_result_mat_1", "five_runif_means_25")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p2 <- round(1-MeanOutFun("five_runif_result_mat_2", "five_runif_means_50")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p3 <- round(1-MeanOutFun("five_runif_result_mat_3", "five_runif_means_75")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p4 <- round(1-MeanOutFun("five_runif_result_mat_4", "five_runif_means_10")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p5 <- round(1-MeanOutFun("five_runif_result_mat_5", "five_runif_means_125")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)

n <- t(rbind(n1, n2, n3, n4, n5))
p <- t(rbind(p1, p2, p3, p4, p5))

tab.2.b <- rbind(n[1,], p[1,], n[2,], p[2,], n[3,], p[3,], n[4,], p[4,], n[5,], p[5,], n[6,], p[6,])

#
p <- 9

n1 <- round(1-MeanOutFun("nine_runif_result_mat_1", "nine_runif_means_25")[c(1:2, 6, 3:5), c(2)]/p, 3)
n2 <- round(1-MeanOutFun("nine_runif_result_mat_2", "nine_runif_means_50")[c(1:2, 6, 3:5), c(2)]/p, 3)
n3 <- round(1-MeanOutFun("nine_runif_result_mat_3", "nine_runif_means_75")[c(1:2, 6, 3:5), c(2)]/p, 3)
n4 <- round(1-MeanOutFun("nine_runif_result_mat_4", "nine_runif_means_10")[c(1:2, 6, 3:5), c(2)]/p, 3)
n5 <- round(1-MeanOutFun("nine_runif_result_mat_5", "nine_runif_means_125")[c(1:2, 6, 3:5), c(2)]/p, 3)

p1 <- round(1-MeanOutFun("nine_runif_result_mat_1", "nine_runif_means_25")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p2 <- round(1-MeanOutFun("nine_runif_result_mat_2", "nine_runif_means_50")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p3 <- round(1-MeanOutFun("nine_runif_result_mat_3", "nine_runif_means_75")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p4 <- round(1-MeanOutFun("nine_runif_result_mat_4", "nine_runif_means_10")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p5 <- round(1-MeanOutFun("nine_runif_result_mat_5", "nine_runif_means_125")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)

n <- t(rbind(n1, n2, n3, n4, n5))
p <- t(rbind(p1, p2, p3, p4, p5))

tab.2.c <- rbind(n[1,], p[1,], n[2,], p[2,], n[3,], p[3,], n[4,], p[4,], n[5,], p[5,], n[6,], p[6,])

table.2 <- rbind(tab.2.a, tab.2.b, tab.2.c)
rownames(table.2) <- rep(rep(c("EBA_Leamer", "EBA_sign", "EBA_90", "MA_unig", "MA_Sala", "MA_BIC"), each=2), 3)
colnames(table.2) <- c(".25", ".50", ".75", "1.00", "1.25")

write.table(table.2, file="pluemper_traunmueller_table2.csv", sep=",")

print(table.2)
##################################################################################################################################
